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ABSTRACT 

A  control  volume  method  is  proposed  for  planar  div-curl  systems.  The  method  is  inde¬ 
pendent  of  potential  and  least  squares  formulations,  and  works  directly  with  the  div-curl 
system.  The  novelty  of  the  technique  lies  in  its  use  of  a  single  local  vector  field  component 
and  two  control  volumes  rather  than  the  other  way  round.  A  discrete  vector  field  theory 
comes  quite  naturally  from  this  idea  and  is  developed  in  the  paper.  Error  estimates  are 
proved  for  the  method,  and  other  ramifications  investigated,  ’  ‘ 


AcCt‘:.IO'i  for 

NTIS 

si 

one 

TALi 

□ 

U'-y.r, 

n 

Jus!;:-; 

j; 

By 

Di  t 

/ 

A 

v<  '«  r»;  * 

v  C.'.dts 

Dist 

A\  . :  v 

i  U  ;  Or 

cul 

1This  research  was  supported  by  the  National  Aeronautics  and  Space  Administration  under  NASA  Con¬ 
tract  No.  NAS1-18605  while  the  author  was  in  residence  at  the  Institute  for  Computer  Applications  in 
Science  and  Engineering  (ICASE),  NASA  Langley  Research  Center,  Hampton,  VA  23665. 

3This  work  was  supported  by  the  United  States  Air  Force  under  grant  AFOSR  89-0359. 


l 


1  Introduction.  Although  div-curl  systems  occur  in  fluid  dynamics,  in  elec¬ 
trodynamics  and  several  other  applications,  relatively  few  discretizations  are 
available.  A  possible  reason  is  that  equivalent  potential  formulations  often  ex¬ 
ist.  Another  reason  is  that  the  simple  minded  Galerkin  finite  element  approach 
is  not  convergent  in  general.  A  least  squares  formulation  is  the  usual  way  to 
get  a  convergent  finite  element  scheme  for  planar  problems.  The  situation  in 
three  dimensions  is  worse.  For  example,  the  vector  potential  approach  can  have 
spurious  mode  problems  [10].  Underlying  the  difficulties  is  the  latent  overdeter¬ 
mination  of  the  div-curl  system  (4  equations,  3  unknown  functions). 

This  article  contains  an  alternative  approach  which  works  directly  with  the  div- 
curl  equations  and  does  not  involve  potentials  or  least  squares.  The  new  recipe 
is  finite  volume  based,  but  differs  in  a  key  ingredient.  The  change  permits 
the  development  of  a  discrete  vector  analysis.  In  turn,  this  provides  tools  for 
analysis  of  the  discretization  and  for  other  purposes. 

Central  to  the  approach  is  the  use  of  dual  pairs  of  meshes  made  up  of  “comple¬ 
mentary  volumes”.  In  three  dimensions,  they  have  the  property  that  the  edges 
of  each  mesh  are  perpendicular  to  the  faces  of  the  other.  There  are  many  possi¬ 
bilities  for  such  mesh  pairs.  In  two  dimensions,  the  simplest  example  consists  of 
the  staggered  Cartesian  meshes  (MAC  meshes)  well  known  in  fluid  mechanics. 
The  complementary  volumes  are  the  basic  mesh  squares  and  the  shifted  mesh 
squares  centered  on  the  nodes  of  the  basic  mesh.  For  triangular  and  tetrahedral 
meshes,  an  example  is  given  by  Voronoi-Delaunay  mesh  pairs.  Many  other  pos¬ 
sibilities  exist,  including  prismatic  meshes  in  three  dimensions  and  combinations 
of  these  meshes. 

The  basic  idea  of  the  discretization  is  to  define  field  components  along  the  edges 
of  one  of  the  meshes,  and,  therefore,  normal  to  the  faces  of  the  other.  In  two 
dimensions,  this  single  component  is  enough  to  permit  the  definition  of  two 
field  operators  corresponding  to  div  and  curl.  With  boundary  conditions,  these 
are  sufficient  to  define  the  discrete  field.  Associated  with  the  nodes  of  the  two 
meshes  are  the  two  discrete  potentials  which  generate  the  null  spaces  of  the 
discrete  div  and  curl.  The  usual  relations  are  valid  between  these  operators  and 
potentials. 

We  will  address  only  the  two  dimensional  problem  in  this  paper.  The  main  ideas 
go  over  naturally  to  three  dimensions  but  are  sufficiently  different  to  warrant  a 
separate  treatment.  It  will  be  given  in  a  forthcoming  report.  The  goals  here  are 
to  provide  a  framework  for  error  estimation  of  complementary  volume  schemes, 
and  to  develop  the  main  tools  of  the  discrete  vector  field  theory. 

In  sections  2  and  3  the  class  of  meshes  of  interest  is  specified,  and  a  typical 
discretization  is  derived.  Sections  4  and  5  contain  the  formulation  and  prop¬ 
erties  of  the  discrete  vector  field  theory  and  some  applications.  Sections  6-8 
are  concerned  with  the  main  results  of  the  error  analysis.  Section  7  in  partic¬ 
ular  establishes  a  link  with  finite  elements  and  applies  the  previous  results  to 
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error  estimate  a  potential  theory  problem.  Section  9  discusses  special  topics 
concerning  rectangular  meshes.  Section  10  extends  the  results  to  more  general 
boundary  conditions. 

The  discretization  technique  reported  here  has  significant  extensions  to  higher 
order  systems  of  partial  differential  equations,  particularly  to  viscous  fluid  flow 
problems.  Algorithms  for  the  Navier-Stokes  equations  are  provided  in  [4]  and 

[7]- 


2  Locally  Equiangular  Meshes.  We  begin  by  defining  a  class  of  meshes  of 
interest.  Let  denote  a  bounded  polygonal  region  of  R? .  Q  may  be  multiply 
connected.  Let  To  denote  the  outer  boundary,  let  i  =  1,2, . .  .,r  denote  the 
inner  (polygonal)  boundaries  if  they  exist  and  let  T  :=  U[_or,.  T;  i  —  1, 2, . . . ,  r 
are  considered  to  bound  “holes”  i  =  1, 2, . . . ,  r  .  Q  will  be  triangulated  and 
a  dual  tesselation  will  also  be  used.  Let  r  denote  a  triangulation  of  Q  with  T 
triangles  denoted  by  r,  i  =  1,2,  ...,T  with  N  nodes  Xi  i  =  1,2, . . . ,  fVj  G  fl 
and  Xj  j  =  Ni  +  l,N\  +  2, . . . ,  N  G  T  as  well  as  E  edges  cr<  i  =  1,2, . . . ,  E\ 
€  fi  ( interior  edges)  and  <x,  =  E\  +  1,  E\  +  2, . . . ,  E  G  T  ( boundary  edges). 

Two  triangles  are  called  adjacent  if  they  share  a  common  side.  Then  we  can 
construct  a  dual  tesselation  by  joining  the  circumcenters  of  adjacent  triangles. 
Many  duals  of  a  given  triangulation  can  be  constructed.  For  example,  another 
one  could  be  made  by  joining  the  centroids  of  adjacent  triangles.  The  dual  which 
is  based  on  circumcenters  is  rather  special  and  will  be  called  the  normal  dual 
since  if  two  triangles  are  adjacent,  the  line  joining  their  circumcenters  is  normal 
to  (and  bisects)  their  common  side.  The  dual  figures  are  polygons  and  in  general 
they  can  have  self  intersecting  boundaries.  They  will  be  called  covolumes.  The 
covolume  associated  with  an  interior  node  is  the  polygonal  figure  obtained  by 
joining  (in  order)  the  circumcenters  of  the  adjacent  triangles  which  share  it. 

We  will  associate  a  (boundary)  covolume  with  each  boundary  node  as  well.  The 
procedure  is  illustrated  in  Figure  1,  where  the  covolume  for  the  boundary  node 
A  is  the  interior  of  the  polygon  PATSRQ. 

In  Figure  1  Q,R,  and  S  are  the  circumcenters  of  their  triangles,  and  P  and  T 
are  the  midpoints  (and  circumcenters)  of  their  edges  wh:  ">  on  T  . 

The  normal  dual,  consisting  of  T  nodes,  E  edges  and  jV  t  v.  i  mes  is  denoted 
by  r'.  Sometimes,  we  will  use  the  “co”  prefix  to  denote  vai.  elements  of  the 
normal  dual,  for  example  in  referring  to  coedges,  comesh  and  so  on. 

There  is  additional  complexity  associated  with  self  intersecting  covolumes  which 
we  wish  to  avoid.  To  do  so,  we  require  that  r  is  “locally  equiangular”  [8]: 

Definition.  A  triangulaiion  is  locally  equiangular  iff  for  every  pair  of  adjacent 
triangles  which  form  a  convex  quadrilateral,  the  sum  of  the  angles  opposite  the 
common  side  is  at  most  180  deg. 
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Elementary  geometry  shows  that  if  r  is  locally  equiangular  then  (1)  each  interior 
covolume  is  convex  and  (2)  distinct  covolumes  have  empty  intersection  and  each 
point  of  Cl  is  either  in  a  covolume  or  on  the  common  boundary  of  two  covolumes. 
Then  r'  has  inner  and  outer  coboundaries  T'  consisting  of  those  coedges  which 
intersect  edges  of  r  with  just  one  node  on  T.  These  will  be  denoted  by  r'0,  and 
T'  i  =  1,2, . . . ,  r  when  the  latter  exist.  The  section  QRS  in  Figure  1  is  part  of 
r'.  Note  that  any  part  of  V  might  penetrate  the  corresponding  part  of  T  .  This 
will  happen  if  there  are  obtuse  angles  opposite  an  edge  of  T,  even  if  r  is  locally 
equiangular. 

We  will  obtain  the  results  assuming  that  r  >  0.  The  modifications  for  r  =  0  are 
mostly  self  evident. 

The  two  tesselations  r  and  r'  are  close  to  being  a  Delaunay- Voronoi  pair.  How¬ 
ever,  although  a  Delaunay  triangulation  is  always  locally  equiangular,  the  con¬ 
verse  is  false,  at  least  when  the  classical  definition  of  the  Voronoi  diagram  is 
used.  The  problem  occurs  at  boundaries.  Recall  that  a  standard  Delaunay  tri- 
angulation  is  defined  by  joining  adjacent  vertices  if  their  Voronoi  figures  share 
a  common  edge.  Then  it  follows  that  the  Delaunay  triangulation  is  a  trian¬ 
gulation  of  the  convex  hull  of  the  vertices.  Since  Cl  is  not  convex  in  general, 
constructing  the  Delaunay  triangulation  of  the  vertices  of  r  does  not  necessarily 
give  a  triangulation  of  Q  .  For  the  purposes  of  this  article  these  distinctions  are 
unimportant.  Local  equiangularity  is  the  only  property  we  need  since  all  the 
properties  required  can  be  derived  from  it. 

Some  of  the  results  below  must  be  modified  for  certain  kinds  of  trivial  triangula¬ 
tions,  in  particular  those  with  no  interior  triangles.  We  will  always  assume  that 
the  triangulations  are  sufficiently  fine  for  the  purposes  at  hand.  No  significant 
loss  of  generality  is  incurred  by  this  assumption. 

We  will  make  frequent  use  of  the  following  special  case  of  Euler’s  formula: 
Lemma  2.1.  Let  r  denote  a  plane  triangulation  with  N  vertices,  T  triangles, 
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E  edges,  and  r  holes.  Then 


N  +T  =  £  +  1  -  r. 


Proof.  If  r  =  0  the  lemma  is  easily  proved  by  deleting  triangles  from  r  while 
maintaining  a  count  of  N,  T  and  E.  If  r  >  0  we  can  imagine  the  holes  trian¬ 
gulated  consistently  with  r  .  This  gives  a  triangulation  say  f  without  holes. 
Combining  the  already  established  result  for  f  and  for  the  triangulations  of  the 
holes,  the  result  follows  by  subtraction.  □ 


3  Div-Curl  Systems.  One  div-curl  system  of  interest  is 


div  u  =  p 

in 

(1) 

curl  u  =  u> 

in  Q 

(2) 

u  -  n  =  / 

on  r 

(3) 

=  7,  i  = 

1,2,. ...r 

(4) 

where  u  :=  («,  v),  curl  u  :=  vx  —  uy,  t  denotes  the  positively  oriented  unit 
tangent,  n  the  outward  unit  normal,  and  it  is  assumed  that 


I  pdxdy  =  I  f  ds. 

J  n  Jr 


(5) 


Also  assumed  is  that  p  £  L2(Q),  u  £  L2(ft)  and  /  £  //1/2(T)  and  that  the 
system  has  a  unique  solution  ’’  £  See  [6], [9]  for  information  on  this 

point.  We  will  explain  the  basic  ideas  of  the  discretization  in  terms  of  (l)-(4). 
Section  10  extends  the  results  to  other  boundary  conditions. 

Referring  to  Figure  2,  we  approximate  (1)  by 

u^hx  +  u2h7  +  u3h3  =  /  pdxdy.  (6) 

Jn 

Here  and  below,  the  Uj  denote  approximations  to  u  ny,  ny  denote  unit  normals 
directed  outwards  and  hj  >  0  denote  the  ordinary  side  lengths  of  r.  There  will 
be  a  similar  equation  for  each  one  of  the  T  triangles  in  r  .  In  matrix  form,  with 
u  denoting  the  vector  of  components  uy,  these  flux  equations  can  be  written  as 

Fu  =  p  (7) 

where  u  £  RE  and  p  £  RT .  Analogous  to  (6)  discrete  fluxes  from  the  holes  are 
also  defined  and  we  denote  them  by 

F(u  Sj)  j  =  1,2, ....  r. 
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For  any  u  G  RE ,  it  is  also  convenient  to  introduce  the  r  x  E  matrix  T  and 
denote  the  hole  fluxes  by  Tu. 

To  approximate  (2)  we  integrate  it  over  an  interior  covolume  r-  £  r'  as  illus¬ 
trated  in  Figure  3.  The  arrows  on  the  covolume  edges  denote  the  directions  of 


the  unit  normals  to  the  associated  triangle  edges.  Approximation  of  the  integral 

I  u  tds  {—l  curl  u  dxdy) 

JdT[  Jr[ 

where  dr^  denotes  the  positively  oriented  boundary  and  t  the  unit  tangent  gives 


Y,Ukh'k  =  f 


u  dx  dy 


(8) 


where  the  sum  is  over  the  coedges  of  dr-  and  h'k  >  0  denotes  the  length  of  a 
coedge.  Assembling  these  circulation  equations  gives  the  matrix  equation 

C0«  =  w  (9) 
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where  u  E  R *  and  u>  E  RNl  ■  The  zero  subscript  is  intended  to  suggest  that 
circulations  in  (9)  are  computed  for  interior  nodes  only.  Discrete  circulations 
around  the  hole  coboundaries  are  defined  similarly  to  (8)  and  denoted  by 

C(u;r')  j  =  1, 2, r. 

These  circulations  will  be  represented  as  Cu  where  u  6  RE  and  C  is  r  x  E.  Ti.e 
boundary  condition  (3)  is  discretized  by  defining  boundary  edge  values  for  u  by 

Ufc  =  f  fds  k  =  El  +  (10) 

hk  Jok 

Note  that  there  are  N  —  N\  of  these  equations. 

The  prescribed  circulations  (4)  are  approximated  as 


9j  j  =  1, 2, ...,  r 


where  sj  denotes  the  strip  lying  between  and  T'- .  We  will  assume  that  T'-  does 
not  penetrate  to  avoid  extending  u  outside  Q  .  This  will  hold  iff  opposite 
every  triangle  edge  on  Tj  j  =  1,2, . . .  ,r  is  an  acute  angle.  (No  restriction  is 
necessary  for  a  simply  connected  domain).  There  is,  of  course,  no  assurance  that 
a  locally  equiangular  or  even  a  Delaunay  triangulation  has  no  obtuse  angles. 

Equations  (7), (9), (10)  and  (11)  are  a  linear  system  of  T  +  N\  +  (N  -  N\)  +  r  = 
T  +  +  r  equations  in  E  unknowns.  By  Euler’s  formula  there  is  one  more 

equation  than  unknowns,  so  we  expect  a  single  constraint  on  the  data.  This  will 
turn  out  to  be 

X>=£/>  (12) 

t.€t  <7,er 

which  holds  by  (5).  We  will  prove  in  section  5  that  these  discrete  equations 
indeed  have  a  unique  solution. 


4  Mesh  Matrices.  In  this  section  some  basic  properties  of  dual  mesh  systems 
are  derived.  These  will  lead  in  the  next  section  to  a  more  detailed  formulation 
of  the  equations  in  section  3.  Similar  results  are  valid  for  more  general  mesh 
systems  than  plane  triangulations  but  will  not  be  needed  below. 

Let  t  denote  an  arbitrary  triangulation  of  D  with  T  triangles,  E  edges  including 
E\  interior  edges  and  N  nodes  of  which  N\  are  interior  nodes.  Label  the  interior 
nodes  1,2,...,^,  the  interior  edges  1, 2, . . . ,  E\  and  assign  the  positive  direc¬ 
tion  along  each  edge  <7*  k  =  1, 2, . . . ,  E  to  be  from  lower  to  higher  node  number. 
Let  r'  denote  a  dual  mesh,  with  T  nodes  E  edges  and  N  covolumes.  The  dual 
can  be  quite  general,  subject  to  having  exactly  one  node  in  each  triangle.  The 
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dual  edges  obtained  by  joining  adjacent  dual  nodes  are  in  a  biunique  correspo- 
nence  with  the  edges  of  r.  For  each  boundary  triangle,  we  join  the  dual  node 
to  a  point  on  the  triangle’s  boundary.  The  dual  edge  corresponding  to  is 
denoted  by  a'k.  For  orienting  the  edges  of  t'  we  use  the  convention  thet  {cr'k.<Tk) 
are  oriented  like  a  Cartesian  coordinate  system  in  the  plane.  This  convention 
applies  to  both  interior  and  boundary  edges. 

Denote  by  D  the  £xiV  edge-node  incidence  matrix  of  r,  where 

{+1  if  is  directed  into  node  j 
—  1  if  <T\  is  directed  out  of  node  j 
0  <Ji  does  not  meet  node  j. 

D  is  easily  seen  to  have  rank  N  —  1.  Define  Do  to  be  the  E\  x  N\  matrix 
obtained  by  deleting  rows  and  columns  of  D  corresponding  to  boundary  edges 
and  boundary  nodes  respectively.  Do  has  rank  N\.  If  r  >  0  we  also  define  an 
Ei  x  r  matrix  V  as  follows. 

{+1  if  (j,  is  directed  into  T5 
—  1  if  is  directed  out  of  T, 

0  c i  does  not  meet  T,. 


In  this,  <r,  denotes  interior  edges,  and  s  =  1,2, . . . ,  r .  V  has  rank  r . 
For  the  dual  r'  we  define  the  E  x  T  incidence  matrix  B  as  follows: 


Bij 


■rl  if  a-  is  directed  out  of  ry 
—  1  if  <r<  is  directed  into  ry 
0  <r,  does  not  meet  r; . 


B  has  rank  T.  Bq  denotes  B  with  rows  corresponding  to  boundary  edges  deleted. 
Bq  is  of  order  Ei  x  T.  If  r  >  0  we  define  an  E  x  r  matrix  B  by 


Bit  := 


+  1  if  cr-  is  directed  out  of  fi, 
—  1  if  <t'-  is  directed  into  fi, 

0  a\  does  not  meet  T, . 


B  has  rank  r. 

An  important  result  is  the  following: 


Theorem  4.1.  Lei  v  G  REl  •  Then  3  <t>  £  RT  such  that  v  ~  Bq4>  tff  D*0v  =  0 
and  Vtv  —  0. 

Proof.  Since  rank  (So)  =  T  —  1,  we  have 

dim  N(Bo)  ~ 
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Ei  -  T  +  1 

Ni+r 


using  lemma  2.1.  Here  and  below,  N(-)  and  /?(•)  denote  the  null  space  and  range 
of  their  arguments.  Next,  it  can  be  verified  directly  that  BqD0iI>  =  0  V  ^  G  Rn'  , 
and  BqT>£  =  0V£  G  Rr ,  so  that  R(D0)  C  JV(B£)  and  R(V)  C  ) .  Direct 
verification  also  shows  that  R(D0)  n  R(V)  =  0.  Since  dim  R(Do)  -f  dim  R( V)  = 
Afj  +  r  it  follows  that  N(Bq)  =  R(Dq)  U  R{T>).  Solvability  of  the  equation 

B0<p  -  v 

holds  iff 

(v,z)=0  Vz  G  N(Bq), 

(  ,  )  denoting  the  standard  Euclidean  inner  product.  By  the  above,  this  is 
equivalent  to 

( v,D0xl> )  =  0  V  V’  6  RN' 

(v,V{)  =  0  G  Rr , 

and  these  in  turn  are  equivalent  to  the  theorem.  □ 

Corresponding  to  theorem  4.1  is: 

Theorem  4.2.  Let  v  G  RE ■  Then  3rl>  £  RN  such  that  v  =  DxL  iff  B*v  0  and 
B‘v  =  0. 

Proof.  Similar  to  theorem  4.1.  □ 


5  Discrete  vector  fields.  In  this  section  we  will  express  the  flux  and  circu¬ 
lation  operators  in  terms  of  the  mesh  properties  of  section  4  and  develop  some 
analogs  of  vector  field  theorems.  For  the  remainder  of  the  paper  we  will  be 
using  the  normal  dual  exclusively.  In  any  normal  dual,  the  distance  between 
two  circumcenters  will  be  zero  if  they  coincide.  Although  this  situation  is  “non¬ 
generic”,  it  does  occur  in  some  situations  (section  9).  Other  than  this,  it  is 
assumed  throughout  that  the  circumcenters  are  all  distinct.  If,  in  fact,  there 
are  coincident  circumcenters,  the  results  obtained  remain  true,  but  minor  varia¬ 
tions  in  some  proofs  may  be  required.  A  second  assumption,  made  throughout, 
is  that  the  coboundary  T'  does  not  penetrate  T.  This  was  already  mentioned  in 
section  3  following  (11). 

In  Re  introduce  the  inner  product  [•,  •]  defined  by 

(13) 


In  (13)  H  :=  diag(/i*)  and  H1  :=  diag(/i't).  Both  H  and  H'  are  invertible  and 
we  let  W  :=  HH'.  In  Figure  4,  ABC  and  ACD  are  adjacent  triangles  from  r, 


[u,  t>]  :=  (u,  HH'v)  =  (u,  H' Hv) 
and  denote  the  associated  norm  by 

IMk=[«,  u]1/2. 
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Figure  4: 


and  P  and  Q  are  their  circumcenters.  Let  |  AC  |:=  />*  and  |  PQ  |:=  h'k .  The 
area  of  the  kite  shaped  figure  APCQ  is  h^h^/2.  The  corresponding  areas  at 
boundaries  are  h\h'J4.  It  follows  that  ||  ||w  is  twice  a  discrete  Lr{Cl)  norm. 

This  interpretation  holds  for  any  locally  equiangular  triangulation. 

Denote  RE  equipped  with  {•,  •]  by  U  ~  Then  we  can  refer  to  “grid 

functions”  u  £  C(Q)  and  regard  them  as  having  boundary  values  u|p  and 
interior  values  u^-  Define 

C'o  '■=  {«  €  U  ;  u|r  =  0}. 

The  flux  and  circulation  operators  of  section  3  can  be  expressed  as 

F  =  BXH 
Co  =  D'0H'0 
C.  =  V'H'0 
T  =  B‘ II 

where  Hg  denotes  the  restriction  of  H'  to  interior  edges..  Verification  of  these 
is  by  direct  calculation.  Note  that  orientations  are  automatically  taken  into 
account  in  this  formulation.  Difference  operators  are  defined  as  follows: 


s<t> 

:=  (//')-'£?<> 

v<>e  RT 

So<h 

:=  (H'0)~1Bo4> 

v<!>£  Rt 

Til ■ 

:=  If~l  Dil’ 

hi’  £  Rn. 

We  will  the  notation  H a  to  denote  the  restriction  of  If  to  interior  edges,  and 
Wo  H0H'o  =  HqHo-  The  theorems  of  section  4  translate  into  theorems  about 
the  existence  of  potentials  for  mesh  functions  u  with  Cou  =  0  and  Cu  =  0 
(velocity/scalar  potential)  and  Fu  =  0  and  IFu  =  0  (stream  function/vector 
potential). 

Theorem  5.1.  If  CqU  ~  0  and  Cu  =  0  then  there  exists  d  £  RT  such  that 
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u  =  5o0.  Conversely,  if  u  :=  5o0  then  Cqu  =  0  and  Cu  =  0. 

Proof.  Use  theorem  4.1  to  prove  the  first  part  and  direct  substitution  for  the 
converse.  □ 

For  the  stream  function,  the  analogous  result  is: 

Theorem  5.2.  If  Fu  —  0  and  Tu  =  0  then  there  exists  if  £  R‘W  such  that 
u  =  Tip.  Conversely,  if  u  :=  Tif  then  Fu  =  0  and  Tu  =  0. 

Proof.  For  the  first  part  use  theorem  4.2.  The  converse  is  easily  proved  by 
direct  substitution.  □ 

Note  that  the  converses  in  these  theorems  furnish  analogs  to  the  vector  identities 
“curl  grad  4>  =  0”  and  “divcurlu  =  0”. 

Another  useful  result  is: 

Lemma  5.1.  For  all  u  G  U ,  <j>  £  RT  and  if  G  RN  we  have 

(1)  [u,  50]  =  (Fu,  4>) 

(2)  [u,  Tif)  =  (Cu,  if). 


Proof.  (1).  By  definition  of  5 

[u,  50]  =  (u,  HH'{H')'xB<f) 

=  (B*H u,  <f)  =  (Fu,  <f). 

The  proof  of  (2)  is  similar. □ 

These  are  analogs  of  integration  formulas.  For  example,  the  first  is  analogous 
to 

/  u.  sj<pdxdy=  —  /  (pdivudxdy  (<f  |r=  0). 

J  n  J  n 

Another  useful  result  following  from  above  is  an  analog  of  the  Helmholtz  de¬ 
composition  of  vector  fields.  In  general,  this  decomposition  is  for  the  subspace 
of  Uq  C  Uq  whose  definition  is: 

Co  :=  {ti  G  Uo,Cu  =  0  and  Tu  —  0} 

Actually,  since  u  G  Uo  =>  Tu  =  0,  this  requirement  in  the  definition  is  redun¬ 
dant.  It  is  included  to  make  explicit  a  symmetry  in  the  hypotheses.  Also,  before 
computing  Cu  we  must  restrict  u  to  interior  edges.  This  point  will  arise  several 
times  below  and  also  in  connection  with  Co-  We  will  still  denote  the  restriction 
by  u,  and  without  explicit  mention  each  time.  If  r  =  0,  then  U0  =  Uo  and  the 
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decomposition  is  for  Uo  itself.  We  also  define  the  following  subspaces  of  rU: 

Z0  =  €  Uo',  Fu  =  0} 

Wo  —  {u  €  Cfoi  Cou  —  0}- 

Now  we  have: 

Theorem  5.3.  Uo  has  the  decomposition 

Uo  —  Zo  ®  Wo, 

which,  is  orthogonal  relative  to  [•,  •]. 

Proof.  First,  note  that  if  to  €  Wo,  then  by  theorem  5.1  and  part  (1)  of  lemma 
5.1VzeUo, 

[*,  tu]  =  [z,  So<t>)  =  [*,  S4>]  =  (Fz,  <6)  e  Rt. 

Hence,  if  z  is  orthogonal  to  Wo,  then  Fz  —  0  which  implies  z  €  Zo.  On  the 
other  hand,  if  some  z  €  Uo  satisfies  Fz  =  0  and  also  Coz  =  0  then  we  have  for 
some  < i>'  €  RT 

[z,  z]  =  [z,  Sd>'] 

=  (Fz,#) 

=  0 


so  that  z  —  0.  This  proves  the  result.  O 

As  an  application  of  theorem  5.4  we  will  prove  that  the  discrete  system  of  section 
3  has  a  unique  solution.  These  equations  are 


Fu  ~  Bx  Hu  =  p 

(14) 

C0u  =  D'qH'qU  ~  u> 

(15) 

Cu~VH'0u  =  g 

(16) 

lu  =  /. 

(17) 

Here,  u  G  RE,  p,  w,  f  and  g  are  as  before  and  the  E  x  (N  —  Ni)  matrix  I 
denotes  the  identity  restricted  to  the  boundary  edge  values. 

Theorem  5.4.  Equations  (14)-(17)  have  a  unique  solution  u  £  U. 

Proof.  Consider  the  homogeneous  problem  Fu  =  0  Cqu  —  0  Cu  =  0  2u  =  0. 
In  this  case,  u  £  To  (1  Zo  and  by  theorem  5.3  it  follows  that  u  =  0,  giving 
uniqueness.  Consequently,  the  matrix  (F(  Cq  I‘  C*]  which  is  of  order  (£7+  1)  x  E 
has  rank  E.  The  augmented  matrix 

'  F ‘  Cf,  I *  O  V 

.  P*  w*  /'  g1  . 
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of  order  (E  +  1)  x  (E  +  1)  has  rank  at  most  E.  This  follows  by  subtraction  of 
the  third  block  row  from  the  sum  of  the  first  block  of  rows  followed  by  use  of 
the  compatibility  condition  (3.12).  Hence,  existence  follows.  □. 

A  method  for  solving  ( 14)-(  17)  by  reducing  them  to  potential  equations  can 
also  be  found  using  the  results  above  including  theorem  5.3.  We  can  suppose 
without  loss  of  generality  that  /  =  0,  simply  by  substituting  the  known  values 
of  u  into  (14)  and  (15).  These  values  do  not  appear  in  (16).  We  can  also  assume 
g  =  0  since  it  is  trivial  to  find  a  solution  of  (16)  and  then  modify  u  accordingly. 
Hence,  it  is  only  necessary  to  consider  (14)-(17)  with  /  =  0  and  g  —  0.  That 
means  u  G  Uo ■  Now  seek  u  as  a  sum  where  both  of  these  are  in  Uo 

and 

Fu(,)  =  p  (18) 

C0u(,)  =  0 

and 

Fu(s)  =  0 
C0u(j)  =  w. 

Consider  the  first  set  of  equations.  By  theorem  5.1  we  can  write  u ^  =  So<j>  and 
(18)  then  becomes 

fs04>  -  p. 

Taking  account  of  the  boundary  conditions,  this  can  be  factored  into 

[SX1/2][W01/2So]<A  =  P- 


The  coefficient  matrix  is  positive  semidefinite  since  B0  has  rank  T  -  1.  Clearly, 
this  reduction  parallels  a  procedure  for  the  continuous  problem.  The  second 
set  of  equations  can  be  solved  by  a  similar  approach.  The  details  are  given 
(for  a  different  context)  in  section  7.  These  procedures  reduce  the  equations 
( 14)-(  17)  to  the  solution  of  positive  semidefinite  equations.  Various  iterative 
procedures  for  the  solution  of  ( 14)-(  17)  turn  out  to  be  disguised  iterations  for 
these  equations. 

6  Error  analysis.  In  this  section  we  will  estimate  the  error  in  approximating 
the  solution  u  of  the  div-curl  system  (l)-(4)  by  the  solution  u  of  the  covolume 
approximation  ( 14)-(  17). 

To  begin,  we  need  the  following  result: 

Lemma  6.1.  Assume  that  v  G  Wo  and  that  u  G  Uo  and  Fu  =  0.  Then 
[u,  v]  =  0. 

Proof.  In  this  lemma  we  do  not  require  that  u  G  Uo-  By  theorem  4.1, and  using 
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ti|r  =  0  we  have  [u,  u]  -  [u,  So4>]  =  [«,  S<j>]  =  ( Fu ,  <fr)  =  0  using  lemma  5.1(1). 

□ 

Next,  we  introduce  the  following  “mesh  functions”  u ^  £  U  i  =  1,2  computed 
from  the  exact  solution  u  of  the  div-curl  system  as  follows: 

ui1*  :=  /  u  n ds  k=l,2,...,E  (19) 

hk  Jak 

where  ak  is  traversed  positively,  and  n  points  to  the  right  of  crk,  and 

«i2)  ~  TT  j  u  n ds  k=l,2,...,Ei 

nk  J°'k 

u[2)  :=  u[l)  k  =  E(l) +  \,...,E 

where  a'k  is  traversed  positively,  and  n  points  along  a'k. 

Now  denoting  :=  u  —  x/'\  i  —  1,2  we  have 

Fc'1*  =  0  ft1)  6  U0 
C0c(2)  =  0  £<2>  6  U0 

Then  by  lemma  6.1 

[e^,  e<2>]  =  0 

and  consequently,  defining  u  :=  (u(I)  +  u(2>)/2,  it  follows  that 

[u-u,  ti-u]  =  [^(£(1)  +  f(2)),  i(£(1)  +  f(2))] 

=  i [«•  ( 1 )  —  f(2),  ft1)  —  f(2)] 

=  -  «(2),  «cd  _  „ca)]. 

The  right  side  is  independent  of  u,  while  depending  only  on  u  so  this  gives  a 
preliminary  evaluation  of  the  error.  We  can  estimate  the  left  side  in  terms  of 
for  example,  by  means  of 

II  u  —  u  ||tv  =  II  u  -  u(I)  -  (u(2)  -  u(1))/2  ||iv 

>  II  U  -  u(1)  \\w  II  u(2)  -  u(I)  ||iv 

so  that  the  estimate  becomes 


II  U  -  u(1)  \\w  <  II  W<2)  -  U(1)  11^  •  (20) 

Evaluation  of  the  right  side  of  this  proceeds  as  follows:  referring  to  Figure  5  let 
K  denote  a  kite  shaped  domain  with  perpendicular  diagonals  s  and  s'. 
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Figure  5: 


We  shall  assume  the  diagonals  are  aligned  with  a  Cartesian  coordinate  system 
intersecting  at  the  origin  O.  Let  the  diagonal  s  on  the  x-axis  extend  from  x  =  —L 
to  x  =  L  and  let  the  diagonal  s'  on  the  y-axis  extend  from  y  =  —Mi  to  y  =  M2. 
Define  linear  functionals  *  =  1,2  on  H  1(K)  by 


P(1)(v) 

m(2)M 


_i _ 

Mi  +  M2 


vi  (0,y)  dy 


where  v  :=  (vlt  v2)  €  H  1(K).  By  the  trace  theorem,  fi^'X')  !  =  1.2  are 
bounded  on  H^/f)  and  so  is  M^X)  —  M^X’)-  Also,  if  v  is  constant  in  K,  then 
(/ / ^  —  /i^)(v)  =  0.  It  follows  that  there  exists  a  constant  C{K)  such  that 

|  (p(1)  -  y(2))(v)  |<  C(A')  |  v  |i,k 


where  |  •  |i  jf  denotes  the  seminorm. 

A  standard  scale  change  argument  shows  that  C(K)  depends  on  the  ratio 


where  M  :=|  M 1  +  M2  |.  In  terms  of  mesh  geometry,  for  the  kite  associated 
with  cf jt  this  becomes 

max 

Substitution  into  the  right  side  of  the  error  estimate  now  gives 


II  «(2)  -  «(1)  \\w  -  X>(1)-a*(2>)(v)  i2i«i 

k 

<  5ZCIuli,Kl,max(/,fc.  (A*)2) 

it 

<  Cmax(/i2,  h'2)  |  u  |2n 
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where  h  :=  max*  /i*  and  h'  :=  max*  h'k.  This  proves: 

Theorem  6.1.  Assume  that  u  £  Then  with  u  denoting  the  approximate 

solution  and  u ^  computed  from  (19),  we  have  the  estimate 

||  u  -  uw  ||  w  <  C  max(h,  ti)  |  u  |i,n  • 

Note  that  there  is  implicit  dependence  on  the  angles  of  the  triangulation  through 
the  appearance  of  the  dual  mesh  parameter  h' . 

This  estimate  can  be  improved  for  certain  regular  meshes.  The  key  observa¬ 
tion  is  that  the  functional  )  —  p^2\)  vanishes  on  linear  functions  in  K 
if  the  vertical  diagonal  in  Figure  5  is  also  bisected.  This  will  occur  in  a  gen¬ 
eral  triangulation  iff  all  the  triangles  have  the  same  circumradius.  A  complete 
characterization  of  such  triangulations  is  not  known.  It  certainly  includes  the 
standard  uniform  triangulation  of  the  unit  square,  and  the  more  symmetric 
triangulation  in  which  each  mesh  square  of  the  uniform  rectangular  mesh  is 
subdivided  by  its  diagonals.  It  also  includes  triangulations  made  of  equilateral 
triangles.  (The  first  two  of  these  examples  are  brought  within  the  covolume 
framework  in  section  10.)  To  exploit  the  mesh  regularity,  we  assume  now  that 
v  £  H2(/y).  Then  it  follows  that 

|  -  p(2))(v)  |<  C(K)  |  v  \2ik 

and  the  constant  depends  on 


•  max(]gr> 

Proceeding  as  above,  we  obtain  for  the  regular  meshes  the  estimate 

||  u  -  u(I)  ||iy  <  Cmax(h2,  (h'f)  |  u  |2,n  . 

From  the  finite  difference  viewpoint,  this  means  that  the  covolume  scheme  is 
second  order  accurate.  The  analysis  strongly  suggests  that  rapid  changes  in  the 
circumradii  of  adjacent  triangles  should  be  avoided  in  practice.  Similarly,  we 
expect  that  meshes  which  vary  smoothly  in  this  sense  will  yield  better  accuracy. 
)  It  would  be  possible  to  estimate  the  deviation  from  second  order  accuracy  in 

terms  of  the  change  in  adjacent  circumradii,  but  we  omit  this. 


7  Covolume  solution  of  Poisson’s  equation.  An  interesting  application 
of  the  stream  function  can  be  given.  Consider  the  equation 


-A*  =  cj  €  ^2(fi)n^(fi)  (21) 

♦|r  =  0.  (22) 
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Integration  over  a  covolume  Tj  gives 

-AS* 

and  using  the  approximation  7V  for  we  get  after  approximation  of  the  line 
integral 

—C^Tip  =  w. 

In  more  familiar  notation,  referring  to  Figure  6,  a  typical  component  of  this 
equation  is 

<»> 

Boundary  conditions  are  used  in  the  obvious  way.  We  assume  a  locally  equian- 


h, 


Figure  6: 

gular  triangulation.  On  rectangular  meshes,  and  to  a  lesser  extent  on  triangular 
ones,  this  approximation  is  well  known. 

Assembling  the  equations  similar  to  (23),  we  obtain  a  linear  system 

Kip  =  u 

where  K  is  of  order  N\  x  N\.  K  has  the  unexpected  property  of  being  the 
same  as  the  piecewise  linear  finite  element  matrix  for  (21)-(22)  on  r.  That  is, 
if  A,  i  =  1,2, . . . ,  Ni  denote  the  standard  piecewise  linear  basis  functions  for  r, 
then 

Ky  =  /(V^XVA j)dxdy  i,  j  =  1,2, . . . ,  Nx. 

J  n 

A  proof  may  be  found  in  [1].  Another  proof  could  be  based  on  the  potentially 
useful  factoring  of  K  which  follows  by  letting  G  :=  Hq1  Do-  In  terms  of  this  we 
have 

I<  =  G'WqG. 
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This  shows  immediately,  for  example,  that  K  is  positive  definite  since  Do  has 
rank  N 

A  difference  from  the  finite  element  formulation  is  apparent  in  the  formation  of 
the  data  vector  ui.  The  components  of  this  vector  are  formed  by  integration  of  u> 
against  the  characteristic  functions  of  the  covolumes.  This  is  in  contrast  to  the 
integration  against  the  A,  called  for  in  exact  finite  element  theory.  It  is  possible 
that  an  error  estimate  for  the  covolume  scheme  could  be  derived  from  the  finite 
element  estimate  by  the  usual  techniques  for  handling  quadrature  formulas.  On 
the  other  hand  it  seems  quite  natural  to  obtain  an  estimate  within  the  covolume 
framework.  This  can  be  done  as  follows.  Define  u  :=  T'i! .  Then  we  have 

Fu  =  0 
CqU  =  w 
u  |p  —  0 

Cu  =  7i  i  =  1,2,  ...,r 

where  7 ,■  is  computed  as  indicated.  The  exact  solution  of  the  analogous  div-curl 
system  is  u  =  (ui,U2)  where  ui  =  dy'il  and  Based  on  u  define  u f1) 

bv  (19).  Since  td1)  |r=  0  and  =  0  by  theorem  4.2  there  exists  £  RN 

such  that  =  T4>^\  Next,  we  have 

Lemma  7.1. 

xl>w  |<=  *(xi)  t=  1,2,..., AT. 


so  that  (T»)  |*=  |t.  Recalling  that  rank(T)=rank(D)= A  —  1,  it  follows 

that  and  xp(x *)  differ  by  a  constant  which  can  take  to  be  zero.  □ 

Define 

a($,  tf)  :=  /  (S7$)(v'Z)dxdy 

Jn 

The  main  result  is: 
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Theorem  7.1.  Let  i p  denote  the  piecewise  linear  inierpolant  of  ip  ike  covolume 
approximation  to  fH.  Then 


II  ^  ^  ilx.n  <  Cmax(A,  ||2i„. 


Proof.  By  theorem  6.1,  we  can  write 

||  u  -  «(1)  |k  =  ||  Tip  -  Tip{l)  ||w 

<  Cmax(/i,  h')  ||  'it  ||j,n 

and  so 

(CoT(ip  -  V-(1)),  V-  -  d>(1))  <  C[max(h,  h‘)  ||  *  ||2,n]2- 


By  the  lemma,  denoting  by  the  piecewise  linear  interpolant  of  ^  on  r, 
a(ip  -  <£,  ip  -  ty)  <  C[max(/i,  h')  ||  $  ||2,n]2 


and 


a(fp  -  fp  -  )  <  C[max(A,  h')  ||  ||2,n]2  +  a(*  -  -  V). 

Using  approximation  theory  it  follows  that 

a(V>  ~  iff ,  ip  -  V)  <  C[max(/i,  h')  ||  ||2,n]2 

and  the  result  follows  from  Poincare-Friedrichs’  inequality.  □ 

Thus,  the  covolume  scheme  and  the  finite  element  scheme  are  of  the  same  order 
in  the  Hl(Q)  norm. 

8  Tangential  Components.  The  discretizations  of  section  3  produce  approx¬ 
imations  to  the  velocity  components  normal  to  the  triangle  edges.  Tangential 
field  components  are  known  in  the  directions  of  the  comesh  edges.  This  is  a 
basic  characteristic  of  the  covolume  method.  But  in  many  applications  it  is  a 
vector  field  that  is  required  and  not  merely  sets  of  components.  An  example 
occurs  in  connection  with  approximating  convection  terms  in  viscous  flow  prob¬ 
lems  [7],  In  this  section  we  give  an  algorithm  for  computing  a  set  of  tangential 
components  from  a  given  set  of  normal  components.  In  this  way,  a  (discrete) 
vector  field  is  obtained.  The  covolume  scheme  itself  is  not  changed.  Instead,  the 
tangential  components  are  found  from  the  covolume  solution  in  a  “postprocess¬ 
ing”  operation.  Analysis  of  the  error  of  the  tangential  approximations  is  also 
presented  in  this  section. 

We  begin  by  noting  that  the  three  normal  field  components  for  a  given  triangle 
are  too  many  to  determine  a  constant  field  in  the  triangle.  The  following  result 
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gives  the  compatibility  condition  for  the  three  components  to  uniquely  determine 
a  constant  field. 

Lemma  8.1.  Given  three  normal  components  t >*,  i  =  1,2,3  on  triangle  sides 
with  lengths  hi  and  corresponding  unit  outward  normals  n,  where  i  =  1,2,3 
there  exists  a  constant  vector  u  such  that 

u  rij  =  Vi  *  =  1,2,3  (24) 


iff  the  flux 


vihi  +  v2h2  +  1*3/13  =  0. 


u  is  unique. 

Proof.  Equation  (24)  is  a  system  of  three  equations  for  two  components  of  u. 
For  nondegenerate  triangles  r  the  coefficient  matrix  A  :=  (rii)j  clearly  has  rank 
2,  so  it  follows  that  dimAf(A()  =  1.  Now  the  equation 


Aw  =  0 


certainly  has  the  solution  w  =  (hi,  h2,  h3),  since 

0=  /  V(1  )dxdy-  /  nds  =  n\hi  +  n2h2  +  n3h3. 

Jt  Jdr 

Hence,  (24)  will  be  solvable  uniquely  iff 

V\hi  +  v2h2  +  1*3/13  =  0 

which  is  what  we  wanted  to  show.  O 

From  this  lemma  it  follows  that  if  Fu  =  0,  then  we  can  construct  a  piecewise 
constant  vector  field  u  on  r  whose  (continuous)  normal  components  are  given 
by  u.  An  easy  way  to  compute  a  tangential  component  along  an  edge  of  the 
triangulation  is  to  average  the  tangential  components  of  the  constant  fields  in 
the  triangles  sharing  the  edge.  The  remaining  problem  is  then  to  construct 
the  piecewise  constant  field  when  the  flux  is  nonzero  and  lemma  81  does  not 
apply.  There  is  no  unique  way  to  do  this.  For  example,  we  can  take  any  pair 
of  normal  components  and  designate  them  as  the  respective  components  of  the 
sought  vector  in  a  triangle.  This  procedure  will  be  used.  Denoting  the  kite  areas 
associated  with  the  triangle  by  K\  >  K2  >  K3,  we  choose  the  corresponding 
components  Ui  and  u2  associated  with  the  largest  and  second  largest  kite  areas 
to  define  the  vector. 

To  compute  a  constant  interpolating  field  u/,  let  u  :=  (uj,  U2)  and  define 
cos 6  :=  rii  ■  112-  Representing  u /  :=  arini  +  02112  gives  the  equations 

U.  r)(:M:)  <“> 
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which  uniquely  determine  oq,  c*2.  It  can  be  checked  from  this  representation 


that 


2  _  _ 
u/  u/  <  — u  ■  u, 
sin  9 


(26) 


a  result  we  will  use  below. 


Interpolation  of  the  u/’s  across  the  common  edge  of  their  triangles  is  also  subject 
to  a  certain  amount  of  arbitrariness.  One  possibility  is  to  use  the  vector  from 
one  of  the  triangles  and  ignore  the  other.  Another  is  to  weight  the  vectors 
equally,  while  a  third  is  to  use  full  linear  interpolation.  In  this  last  case,  it  is 
necessary  to  assign  a  localization  of  the  vectors.  The  most  natural  choice  is 
to  attach  each  u/  to  the  circumcenter  of  its  triangle.  To  justify  this,  we  can 
observe  that  the  two  directions  from  which  u/  is  made  indeed  pass  through  the 
circumcenter.  The  line  joining  the  circumcenters  is  divided  in  a  certain  ratio  by 
the  common  edge  of  the  triangles  which  determines  the  interpolation  weights  in 
the  usual  way.  Note  that  if  one  of  the  triangles  is  obtuse,  this  becomes  a  linear 
extrapolation.  Each  of  these  methods  corresponds  to  a  weighting  of  the  form 
AiU|^  +  in  which  \\  +  A2  =  1.  Some  situations  suggest  other  weightings 

to  use.  At  boundaries  in  particular,  a  (1,0)  weighting  is  required.  In  other  cases 
the  weights  can  vary  from  edge  to  edge. 

Next,  we  will  estimate  the  error  of  this  tangential  approximation  scheme.  The 
interpolated  vector  for  the  kth  edge  has  the  form 

w<*>  :=  A^u^  +  A^V/’) 

where  the  superscripts  (&i),  (A'2)  denote  the  two  triangles  on  each  side  of  the 
kth  edge.  The  tangential  component  of  this  is  just  •  w^.  For  u  E  U ,  let 
Tu  E  U  denote  the  mapping  to  the  tangential  components  generated  in  this 
way. 

A  parameter  which  will  appear  in  the  estimate  is  the  ratio  \J I\[^ / —■  6^k\ 

where  A'^  are  the  kite  areas  for  rj.  We  will  define  6  maxj  (6^).  This 
quantity  will  be  estimated  later.  Also  we  will  define  0  :=  minm(|sin  0(m,|) 
where  corresponds  to  the  angle  in  (25)  above  and  the  minimum  is  over  the 
triangles.  In  addition  to  these  we  will  denote  by  A  the  maximum  absolute  inter¬ 
polation  coefficient:  A  :=  max*(|Ajfc^|,  lA^I).  The  next  result  will  be  required 
below: 


Lemma  8.1.  The  following  bound  on  T  holds: 


Mw  <  ^-Nk 


Vu  G  U. 


Proof.  Let  denote  the  kite  area  on  the  edge  Then  since 
|t(*)-w(*)|*  <  2A2(||k/‘)||=,+||u^)||fa) 
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we  have,  using  (26) 

]T  \(Tu)k\2K(k)  =  ^  |t(t)  •w(t)|2/\(t) 

t  Jc 

k 

Also 

=  [(«<*■  >)2  +  („<*'  >)2]A(i) 

<  62[(u(1*‘,)2Af‘)  +  (4il))2Kf‘)] 

with  a  similar  bound  for  Noting  that  each  product  appears  at  most 

four  times  and  substituting  these  into  the  previous  inequality  gives  the  result. 

□ 

Next,  let  u  G  H'(Q),  and  define  u{1)  by  (19).  In  addition,  let  i;(1'  be  defined 
by 

t/i1 1  :=  7—  f  tu ds  it  =  1,2 . E  (27) 

hk  Jak 

where  ak  is  traversed  positively  and  t  denotes  the  unit  tangent  along  <rk.  Now 
we  have: 


Lemma  8.2.  The  estimate 

||i>ll)  -Tu(1)||tv  <  —  max(A,^')IMIi,n 

holds. 

Proof.  Let  ak  denote  an  edge  of  Tj  G  r.  The  functional 

uj(u)  :=  (4°  -  tk  ■  u(/}) 

is  linear,  bounded  on  H 1  ( r;)  by  the  trace  theorem,  vanishes  on  constant  u,  and 
so  has  the  estimate 

Mu)l  <  C(7;)IuIi,tj- 

By  the  usual  mapping  to  a  reference  element  [3]  it  follows  that 


C(Tj)  < 


c 

0' 


Following  the  notation  of  the  previous  proof,  using  Ai  +  A2  = 
(r(1)  -  Tu(n)t  <  Adi/jt.fu)!  +  |j/fcj(u)| 

c\ 

<  -0-(l«ll.TM  +  |u|ilTkj). 


1  we  have 
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and  consequently  with  a  modified  constant  C, 

'£\(v™-T^)t\iKk  <  (^)2max(tfW)X>|2,rj 

k  j 

which  is  equivalent  to  the  required  result.  □ 

Combining  lemmas  8.1  and  8.2  now  gives: 

Lemma  8.3.  Let  u  £  U  denote  the  covolume  approximation  to  u  £  H’(f2) 
estimated  in  theorem  6.1  and  let  Tu  denote  the  approximate  tangential  compo¬ 
nents.  Then  the  estimate 

||u(1)  -  Tu\\w  <  ^—mai  (/i,  h^ljuUi.n 

holds  where  t/1*  is  defined  by  (27). 

Proof.  It  is  only  necessary  to  observe  that 

Hi/1)  -Tu\\w  <  Hi/1*  -Tv)^\\w  +  \\T(u  -  ti(1))||w 

and  to  use  lemma  8.2  on  the  first  term  and  lemma  8.1  and  theorem  6.1  on  the 
second.  □ 

We  may  now  construct  the  discrete  vector  field  uc„  :=  (u,  Tu)  £  U  x  U  and 
regard  it  as  a  covolume  approximation  to  u  the  solution  of  the  div-curl  system 
(l)-(4).  The  error  estimate  for  ||(u  —  Tu  —  r^))||2v  :=  ||u  —  1  ^ || ^  +  ||Tu  — 

follows  immediately.  Defining  u(1)  :=  (u^\  t/1)  we  summarize  with 

Theorem  8.1.  The  covolume  approximation  u(1*  satisfies 

l|uc  -  u(1)||vv  <  ^-max(ft,  A')||u||i,n 

where  0  denotes  minm  |sin0m|,  6m  denoting  the  angles  of  the  triangles  in  r,  A 
denotes  the  largest  absolute  interpolation  coefficient,  and  6 2  denotes  the  largest 
kite  ratio  for  triangles  in  the  triangulation. 

Proof.  This  is  simply  a  question  of  applying  the  definition  of  the  norm  on  the 
left  and  theorem  6.1  and  lemma  8  3.  □  . 


We  can  make  this  bound  more  explicit  for  acute  angled  triangulations.  In  that 
case,  we  have  A  <  1.  For  any  acute  angled  triangle  tj  £  r,  let  Ap  >  AP  > 
Ap  >  0  denote  the  areas  of  the  three  triangles  formed  when  the  circumcenter 
is  joined  to  the  vertices.  Thus  AP /\rj\  are  the  barycentric  coordinates  of  the 
circumcenter.  Obviously,  AP  >  |ij|/3-  A  calculation  shows  that  A P  <  2A2;). 
We  will  also  use  the  area  ratio  of  the  triangulation  defined  by 


rT 


max 

■j 


hi 


(r,,  tj  adjacent). 
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Note  that  this  ratio  is  formed  for  adjacent  triangles  only. 

Denote  the  kite  areas  of  Tj  by  A'p*  >  >  0  (since  h'k  >  0).  Below, 

we  will  use  the  fact  that  for  acute  triangles  we  always  have  A[j]  <  1<{J\  This 
follows  since  the  maximal  property  of  ensures  that  it  is  always  at  least  as 
large  as  the  '  which  are  not  part  of  I\[J\  But  at  least  one  of  these  two  A ^ 
must  be  at  least  equal  to  A^ ,  so  the  above  property  does  hold.  It  follows  that 

A'i(^)  =  hkh'k/  2 

<  Mir,)  +  rr|r,| 

<  2y42(r;)  +  rr3.4i(r;) 

<  2Kn(r] )  +  rT3  •  2K2(Tj ) 

<  8  rr/\2(r;). 

This  shows  that  we  can  take  8  =  \/SrT .  It  would  of  course  be  possible  to  bound 
rT  in  terms  of  mesh  lengths  and  0  but  we  prefer  not  to  do  this.  Instead,  making 
use  of  this  value  for  6  in  theorem  8.1  now  gives  the  estimate 

||uc„  -  u(1)||iv  <  C ^  max  (h,  /i')||u||1  n- 
This  is  valid  for  acute  triangulations. 

Estimating  A  and  6  is  more  difficult  when  there  are  obtuse  triangles.  The 
problem  is  associated  with  h'k/hk  values  which  approach  zero.  This  will  not  be 
discussed  further  here,  but  section  9  has  some  related  information. 

9  Comments  on  rectangular  meshes.  The  earlier  results  have  been  ob¬ 
tained  under  the  assumption  that  all  h'k  >  0.  This  is  equivalent  to  assuming 
that  circumcenters  of  adjacent  triangles  are  distinct.  Now  we  will  briefly  summa¬ 
rize  the  necessary  changes  for  the  contrary  case.  The  previous  results  continue 
to  hold  with  some  slight  changes. 

Coincident  circumcenters  will  occur  if  the  triangulation  contains  a  cyclic  quadri¬ 
lateral  Q.  The  distance  between  the  circumcenters  of  the  two  triangles  making 
Q  is  h'  :=  0.  There  is  no  change  in  the  definition  of  the  covolumes,  and  the 
discrete  circulation  equations  (9)  continue  to  approximate  the  continuous  ones 
(2).  Note  that  the  normal  component  for  the  diagonal  does  not  appear  in  the 
circulation  equations.  The  flux  equations  are  also  unchanged  by  the  existence 
of  Q,  but  the  diagonal  component  does  appear  in  them.  A  problem  shows  up 
in  the  anal  .  sis,  since  the  kite  area  for  the  diagonal  is  zero  and  consequently 
||  ■  ||w  is  now  merely  a  seminorm.  This  inconvenience  is  easily  removed  by  a 
slight  modification  of  the  equations:  we  simply  add  the  two  flux  equations  asso¬ 
ciated  with  Q.  The  diagonal  component  does  not  appear  in  the  resulting  sum. 
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In  this  way,  the  number  of  equations  is  reduced  by  one  and  so  is  the  number  of 
equations. 

More  directly,  we  can  modify  the  triangulation  by  deleting  the  diagonals  of  all 
cyclic  quadrilaterals.  Circumcenters  are  well  defined  for  each  of  the  resulting 
mesh  figures.  Flux  and  circulation  equations  are  written  in  the  obvious  way. 
Euler’s  formula,  lemma  2.1  is  still  true  and  there  will  be  a  single  consistency 
condition  analogous  to  (5).  The  efTect  of  this  modification  is  to  remove  the  dual 
edges  which  have  h'  =  0.  With  these  gone,  we  can  prove  the  key  theorem  5.1. 
which  requires  the  existence  of  Still  defining  the  kite  areas  by  joining 

circumcenters  of  mesh  figures  to  their  vertices,  the  results  in  sections  6-7  follow 
essentially  verbatim. 

These  observations  apply  particularly  to  cartesian  rectangular  meshes.  Normal 
to  vertical  (horizontal)  mesh  segments  we  have  horizontal  (vertical)  field  compo¬ 
nents.  Flux  equations  are  written  for  mesh  rectangles,  and  circulation  equations 
for  the  dual  mesh  rectangles  around  the  mesh  points.  The  two  potentials  <f>  and 
V  are  defined  respectively  at  the  mesh  cell  circumcenters  and  the  mesh  points 
themselves.  The  mesh  pair  is  the  same  as  the  usual  staggered  mesh  system 
found  in  fluid  flow  discretizations.  Error  estimation  of  this  scheme  is  given  at 
the  end  of  sect  ion  6,  where  second  order  accuracy  is  demonstrated.  The  scheme 
itself  is  advocated  in  [2]. 

There  are  occasions  when  rectangular  and  triangular  meshes  can  be  usefully 
combined.  This  is  especially  true  when  highly  stretched  meshes  must  be  used. 
In  regions  where  stretching  is  necessary,  rectangular  mesh  cells  can  be  used  in 
place  of  triangles.  This  might  be  helpful  in  calculations  involving  boundary  and 
other  kinds  of  layers. 

A  situation  related  to  the  above  occurs  when  some  h'k/hk  is  small,  but  positive. 
Although,  iri  a  sense,  this  is  inconsequential  from  the  analytical  viewpoint,  in 
finite  precision  environments  instabilities  could  result.  The  three  dimensional 
case  is  worse  in  this  respect,  because  there  is  more  opportunity  for  degeneracy. 
In  practice,  it  should  be  satisfactory  to  merely  set  to  zero  those  b!  for  which 
h'k/hic  is  below  some  threshold  and  proceed  as  indicated  above.  This  does 
introduce  some  extra  error,  but  by  a  suitable  choice  of  the  threshold  it  should 
be  within  the  discretization  error. 


10  Mixed  boundary  conditions.  In  this  section,  we  will  extend  the  earlier 
results  to  boundary  conditions  which  are  partly  normal  and  partly  tangential. 
The  technique  will  be  illustrated  for  the  simply  connected  problem.  If  required, 
multiply  connected  cases  can  be  handled  by  direct  extension  of  the  method 
already  used  in  earlier  sections. 

We  consider  in  the  simply  connected  polygonal  domain  Q  the  system 

div  u  =  p  in  12  (28) 


24 


curl  a  =  u) 

in  Q 

(29) 

u  •  n  =  / 

on  rn 

(30) 

II 

on  r, 

(31) 

where  T  =  Tn  U  T(  and  neither  part  has  zero  measure.  For  illustration, we  can 
assume  that  T  is  divided  into  two  continuous  parts  Tn  and  r<  each  with  nonzero 
length.  We  assume  that  p  €  L2(Q  ),  u>  G  L~(Q  )  and  /  £  Hl^{V  )  and  that  the 
system  has  a  unique  solution  u  £  H!(fi  ) 

Discretization  of  (28)-(31)  follows  section  3,  except  that  (10)  is  applied  only 
to  boundary  edges  on  T„.  Here  and  below,  it  is  assumed  that  nodes  of  the 
triangulation  always  separate  F„  and  I\  and  that,  conventionally,  the  separating 
nodes  belong  to  T,.  Boundary  edges  of  r  which  are  in  rn  are  labelled  E\  + 
1, ...,  £2  and  boundary  edges  in  F*  are  labelled  £2+1, ...,  £.  Then  in  place  of 
(10)  we  have 

Uk  =  —  f  fds  k  =  £1  +  1, . . . ,  £2. 

hk  J<jk 

To  complement  the  boundary  conditions,  some  boundary  circulation  equations 
are  used.  The  derivation  of  these  equations  will  be  illustrated  using  Figure  7. 
ABC  is  part  of  ft  and  the  boundary  covolume  r'B  associated  with  the  node  at 


Figure  7: 

B  is  shown.  The  orientation  conventions  are  unchanged.  Define 

■=mLsds 

and 

,,c  :=  m  L3ds' 

where  AB  and  BC  are  the  positive  directions.  The  circulation  equation  for  dr'B 
is  approximated  by 


u'hi  +  9  yABhA 

9t‘b 


B  +  “  I’BchBC  =  J 


u >  dx  dy 


(32) 
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and  a  similar  equation  is  used  for  each  boundary  covolume  associated  with  a 
node  strictly  interior  to  IV  It  is  always  assumed,  without  explicit  mention, 
that  h  is  sufficiently  small  for  T(  to  contain  at  least  one  strictly  interior  node. 
This  assumption  is  essential  for  what  follows.  If  it  is  not  satisfied,  then  we  are 
essentially  in  the  situation  of  the  previous  sections  as  far  as  the  discretization 
is  concerned. 

The  system  of  equations  which  results  may  be  written  as 


e 

II 

(33) 

CbU  =  w  -  g 

(34) 

u|r  =  f 

(35) 

where  g  is  found  from  the  terms  involving  the  tangential  boundary  function 
in  equations  like  (32)  and  Ct,  has  no  explicit  dependence  on  the  lengths  of  the 
boundary  edges. 

The  number  of  nodes  interior  to  T,  is  E  —  E2  and  the  number  of  edges  is  therefore 
E  —  E2  —  1.  In  all,  the  number  of  equations  is 

T  +  Ni  +  (E  -  Ei  -  l)  +  (E2  -  £1)  =  T  +  N-l 

=  E 

by  lemma  2.1.  This  time,  unlike  the  earlier  situation  the  number  of  equations 
and  unknowns  are  equal. 

Now  we  will  prove  that  the  solution  of  the  linear  system  is  unique.  For  this, 
we  will  generalize  the  Helmholtz  type  theorem  5.3  to  cover  the  new  boundary 
conditions.  It  is  necessary  first  to  have  results  analogous  to  theorems  4.1  and 
5.1. 

Recalling  the  definition  of  the  matrix  D  in  section  4,  let  Db  denote  the  restriction 
of  D  to  {x*;  x;  £  U  rj,  and  let  Sj  denote  B  with  columns  E2  +  1, . . . ,  E 
deleted.  Thus,  dim R(Bb)  =  E\  +  E  —  E2  =:  E' . 

Theorem  10.1  Lei  Q  be  simply  connected  and  assume  that  v  £  RE  .  Then 
there  exists  <t>  £  RT  such  that  v  =  Bb<t>  iff  D\v  =  0. 

Proof.  The  proof  of  this  follows  closely  the  proof  of  theorem  4.1  and  we  shall 
omit  most  of  the  details.  We  have  to  solve 


Bb4>  =  v. 


That 


B\Dbi}>  =  0 


follows  from  the  easily  proved  relation  B(  Dip  =  0.  Also,  rank  (Db)  —  N' ,  where 
N'  denotes  the  number  of  interior  nodes  plus  nodes  on  rt.  In  addition, 


dim N(B\)  =  E'  -T. 
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Counting  and  use  of  lemma  2.1  shows  that  N'  —  E'  —  T  so  the  result  follows.  □ 

Now  we  have  that  C(,u  =  0  iff  D\H'u  =  0.  Define  Sh  :=  (//')_1J3j  (denoting  the 
restriction  of  ( H')~l  to  R{Bb)  by  the  same  symbol).  It  follows  that  CbU  =  0  iff 
u  =  St,4>  for  some  <j>  €  RT  ■ 

Next,  define 

U'Q  :=  {u  €  RE-u\r,  =  0} 

W'Q  :=  {u  G  Uo,CbU  =  0} 

2o  :=  {«  €  t/J;Fti  =  0}. 


Theorem  10.2  Assume  that  Cl  is  simply  connected.  Then  U q  has  the  decompo¬ 
sition 

K  =  W'0®  z'o 

which  is  orthogonal  relative  to  [-,  •]. 

Proof,  w  €  W0  implies  [it,  w]  =  [it,Sj<£]  =  [tt,  S<fi]  =  (Fu,  <j> ).  Hence,  [u,  u>]  = 
0  for  all  w  6  Wq  implies  that  Fu  —  0.  The  rest  of  the  argument  is  as  before.  □ 

Existence  and  uniqueness  for  the  discrete  system  follows  directly,  since  if  u  is 
a  solution  of  the  homogeneous  system  we  have  u  6  Wq  D  Z'q  =  0.  Existence 
follows  from  uniqueness  since  the  coefficient  matrix  is  square. 

It  remains  to  estimate  the  error  of  this  approximation.  The  method  is  similar 
to  section  7  with  some  small  modificaions.  We  define  uO  just  as  in  (19).  Cor¬ 
responding  to  t/2)  we  define  a  very  slightly  different  function,  still  denoted  by 
i/2)  as  follows: 

ujj.2)  :=  7^-  /  u  ■  n  ds  k  =  1, 2, . . . ,  E\  k  =  £2  +  1,  •  •  • ,  E 

hk  Jo'k 

u[2)  :=  u[l)  k  =  +  1,...,£2. 


The  edges  crk  k  =  E^  -f  1, ...,  E  are  in  r(  ;  each  of  the  coedges  (r'k  extends  from 
crjc’s  midpoint  to  the  circumcenter  of  the  triangle  containing  cr*,  exemplified  by 
PQ  and  TU  in  Figure  7.  As  before,  define  =  u  —  u(,)  i  =  1,2  .  Then 

Fe^  =  0  e(1)  €  Z'a 
Cjf(2)  =  0  e(2)  e  W. 0 

and  [flD,  e^2^]  =  0  by  theorem  10.2.  Following  the  method  of  section  t>,  we 
obtain,  corresponding  to  (20) 

||u  —  u^[|(v  <  ||u^2'  — 
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The  right  side  can  be  estimated  as  before,  except  that  now  there  will  be  a 
contribution  to  the  norm  from  the  boundary  segments  along  IV  The  earlier 
estimation  technique  works  here  too,  and  is  especially  simple  if  each  coedge 
associated  with  T(  lies  wholly  in  ft.  This  will  certainly  be  the  case,  for  example,  if 
the  boundary  triangles  are  acute  angled,  since  then  the  coedges  will  lie  inside  the 
triangles  themselves.  That  is  the  simplest  case.  We  summarize  this  discussion 
with: 

Theorem  10.3  For  the  equations  (28)-(31)  and  their  discretization  (33)- (35), 
we  have  the  estimate 

iv  <  Cmax(/i,  /i')|u|i  fi. 


So  far,  we  have  considered  the  purely  normal  boundary  condition  and  the  mixed 
boundary  condition.  For  the  purely  tangential  case  the  simplest  approach  seems 
to  be  to  use  tangential  field  components  in  place  of  the  normal  ones  used  in 
sections  2-6.  The  circulation  equation^  are  formed  using  the  triangle  boundaries 
and  the  flux  equations  are  formed  using  the  covolume  boundaries.  It  is  apparent 
that  this  method  gives  a  natural  dual  to  the  earlier  one.  We  will  not  go  into  the 
details  here,  but  an  almost  verbatim  development  of  the  theory  can  be  carried 
out.  This  includes  most  of  the  results  of  sections  2-6.  The  main  difference  is  to 
exchange  the  roles  of  <j>  and  ip.  This  duality  is  strongly  reminiscent  of  complex 
variable  theory.  In  fact,  a  similar  connection  was  developed  in  [5].  The  results 
given  above  provide  useful  tools  for  the  analysis  of  that  scheme. 
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